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Motivated by the large strain shear of loose granular materials we introduced a model which 
consists of consecutive optimization and restructuring steps leading to a self organization of a density 
field. The extensive connections to other models of statistical phyics are discussed. We investigate 
our model on a hierarchical lattice which allows an exact asymptotic renormalization treatment. 
A surprisingly close analogy is observed between the simulation results on the regular and the 
hierarchical lattices. The dynamics is characterized by the breakdown of ergodicity, by unusual 
system size effects in the development of the average density as well as by the age distribution, the 
latter showing multifractal properties. 
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I. INTRODUCTION 

Slow dynamics with no separation of time scales rep- 
resent a major challenge of statistical physics. Experi- 
mental or simulation approaches are extremely difficult, 
so in most cases new ideas and models are needed for the 
understanding of this kind of problems. 

There can be different roots of slow dynamics: Systems 
close to the critical point slow down enormously due to 
the increasing characteristic time. Phase separation is 
often accompanied by a slow coarsening process Q|. In 
glasses the free energy landscape is so complicated and 
structured that the system never finds the global min- 
imum and shows a history dependent behaviour called 
aging ||. Slow dynamics may also occur in intrinsi- 
cally dynamic, driven systems leading to scale free fractal 
structures. The name of self-organized criticality covers 
a whole family of related models ||. 

In this paper we study a model (introduced in 
where the system exhibits a very slow evolution with a 
tendency of getting stuck in metastable states. However 
the model is different from those studied earlier in the 
sense that there is an element of both energy as well as 
entropy barriers being present as a result of the rules of 
evolution. We are able to directly link the slow evolution 
to a break down of ergodicity in the dynamics. This then 
leads to several interesting features of the model such 
as non-trivial system size effects, a multi-fractal "age" 
distribution and a non-trivial temporal evolution. 

Motivated by our study of shearing loose granular ma- 
terials |||, we report in this paper about a new mech- 
anism leading to slow dynamics. In granular materi- 
als displacement occurs in a localized manner, in 'shear 
bands' which are formed along the weakest parts of the 
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samples. During shear, grains can rearrange themselves 
and occasionally strengthen the local structure. In such 
a case, the shear band finds a new configuration which 
avoids this zone. Based on this picture we introduced 
a model where consecutive steps of optimization (find- 
ing the weakest part of the sample) and restructuring 
(random rearrangement of grains) takes place. Assuming 
translational invariance in the shear direction the model 
becomes two-dimensional. We have studied the model 
numerically in detail on regular lattices Q, however, it is 
difficult to go beyond the simple description of numerical 
simulations. It turned out to show unexpected proper- 
ties including extremely slow dynamics and unusual size 
dependence (breakdown of ergodicity), and it provides 
with interesting predictions for the granular system. 



The aim of this paper is to study this same model 
on a hierarchical diamond lattice both numerically and 
analytically, and to compare these results with the sim- 
ulations on the euclidean lattice. We find that despite 
the very different connectivities of these two lattices, the 
qualitative behaviour is much the same; for some proper- 
ties, there is a quantitative matching as well. The recur- 
sive nature of the hierarchical lattice however aids the 
analytical treatment greatly, thus helping us getting a 
deeper understanding of the problem. 



The paper is organized as follows: In the next section 
we define the model in general and on the hierarchical 
lattice. In Section III, the relation of the model to other 
problems of statistical physics is discussed. In Section 
IV, the numerical results are shown and compared to the 
regular lattice simulations. In Section [y] we present the 
exact asymptotic solution of the model. We conclude in 
Appendix A and B contain technical details 



section VI. 

of the calculations used in Section [y| 
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II. THE MODEL AND THE HIERARCHICAL 
LATTICE 

The model that we study in this paper is defined as fol- 
lows: A two-dimensional field is characterized by a single 
scalar parameter, the density g(x,y). Initially this den- 
sity is generated randomly from the distribution pi(g). 
At every step we search for the minimal path V* that 
is defined as follows: The minimal path is a continuous, 
directed path V that spans the system in the x direction 
and the sum S of the local densities along it, 

S{V)= e(x,y) (i) 

(x,y)£V 

is minimal among all possible paths. The minimal path 
is the path V* for which S(V*) is minimum. 

Once the minimal path is found the density values of 
the points belonging to the minimal path (a;, y) £ V* 
are replaced by new densities randomly picked from the 
distribution p r {g)- 

The above process is repeated as long as desired. A 
single time step consists of both searching for the minimal 
path, as well as refreshing the local densities along it. 

In the following we restrict ourself to the case where 
Pi and p r are uniform distributions in the interval [0 : 1] . 
Our model is discretized on a lattice. In || we report 
on detailed numerical results for various properties of 
the model on the Euclidean square lattice. The analytic 
treatment on the square lattice has not been possible so 
far. However, in this paper, we obtain exact asymptotic 
solutions for the model on the hierarchical diamond lat- 
tice. 

The hierarchical diamond lattice JjJ is constructed as 
follows: We first consider a single bond connecting two 
points A and B. This constitutes the most elementary 
(generation 0) lattice. The first generation lattice is ob- 
tained by substituting the unique bond by an elementary 
"diamond" of four bonds, i.e. two parallel connections 
each consisting of two bonds in series [Fig. [l] b)]. The 
next generation is obtained recursively by the substitu- 
tion of each bond by a diamond [Fig. |l| c)] . Repeating the 
above procedure N times, produces a iVth-generation hi- 
erarchical lattice. This lattice has a dimension equal to 2, 
and hence can be compared to its Euclidian counterpart. 

All results are based on the exploitation of the recursiv- 
ity of the construction of the lattice. If one can compute 
the properties of an elementary diamond and transform 
this into a single bond endowed with the same, a recur- 
sive use of this procedure clearly allows the reduction of 
the entire lattice back to a single bond thus determining 
the global behavior. This is a real space renormalization 
procedure and the structure of the lattice makes such 
renormalization treatments exact. Hierarchical lattices 
have been widely used to study several phenomena such 
as percolation [|| , spin models Q , sums of directed paths 
p0| etc. However, usually there is a price to pay in that 
the result may differ from its Euclidean lattice counter- 
part. There is no general formalism by which means to 
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FIG. 1: The Oth a), 1st b) and 2nd c) generations of the 
hierarchical lattice. 

estimate the validity of hierarchical lattice results for the 
Euclidean lattice. Therefore, it is necessary to resort to 
numerical results to assess the similarity between the two 
cases. It will be shown in the following that indeed the 
analogy between the results obtained on both kinds of 
lattices is extremely close. Therefore, the analytical so- 
lution obtained here provides a better understanding of 
the Euclidean lattice case. 



III. RELATION TO OTHER STATISTICAL 
PHYSICS MODELS 

Before reporting the result of numerical simulations 
on the hierarchical lattice, we point out some analogies 
which can be drawn between our model and other diverse 
problems of statistical physics. 

The rules of our model, finding the extremal directed 
spanning path at every instant, is similar to finding the 
ground state of a directed polymer in a random potential 
jnj. However, in our case this potential is uncorrelated 
only at the beginning; it changes in time through the 
process already described above, of ascribing new den- 
sities to all sites along the minimal path. The shape of 
the path on Euclidean lattices is found to be self-affinc 
in the directed polymer problems. The model studied 
here changes the underlying potential landscape in a self- 
organized way and naturally induces correlations. These 
in turn change the self-affine exponent of the path. This 
feature is studied numerically in an earlier paper j5| . 

The rules of our model can be regarded as a generaliza- 
tion of the Bak-Sneppen model of evolution |12f , but in 
higher dimensions. Indeed, the constraint of finding the 
minimal path and then changing it, puts this model in a 
class of extremal models studied in contexts as different 
as interface depinning |D| and flux creep However 
there is an important difference between our model and 
other extremal models. In the latter case, the system 
(usually one-dimensional) reaches a steady state which is 
'self-organized critical', in the sense that there is a power- 
law distribution for avalanches in the steady state. In 
the case of our model, no steady state is reached and all 
quantities depend on time. As we will see, we can define 
avalanches which are indeed power-law distributed, but 
always with time dependent prefactors. The difference 
is best illustrated if we look at the simplest of these ex- 
tremal models, the Bak-Sneppen model jl2). This model 
is defined on a one dimensional periodic array of ran- 
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dom numbers where at every update, the least and its 
neighbours are refreshed from a given time independent 
probability distribution. Our model is however related to 
a variant of this in which only the least is changed. In one 
dimension, changing only the least does not lead to a very 
interesting behaviour. However in two dimensions, as we 
will see, changing only the minimal path leads to very 
non-trivial effects. Further, the simple minded variation 
of the original Bak-Sneppen model turns out to be very 
useful in solving our model on the hierarchical lattice. 

There are also connections between this model and the 
apparently unrelated problem of a random walk in a dis- 
ordered potential. If we consider a one-dimensional cross- 
section of the model perpendicular to the minimal path, 
we could imagine the point through which the interface 
passes through, as the position of a random walker. The 
subsequent dynamics can then be interpreted as that of 
a walker moving through an initially random potential, 
modifying it along the way. While the actual dynamics 
of the interface in two-dimensions is quite complicated 
to translate in its entirety into one-dimension, it is pos- 
sible to do so in the simplest case when only corner flips 
are allowed for the interface. For the walker, this simply 
translates to the condition that the subsequent position 
of the walker is on one of the neighbouring sites of the 
present one, chosen by an inequality condition. The value 
at the site the walker has just left, is also changed. We 
have studied such an active walker model in detail Q 
and find that it leads to logarithmically slow dynamics. 

Because of the extremal condition used in finding the 
minimal path at each time step, the solution of the model 
on the hierarchical lattice uses results from extreme- value 
statistics jTfJ. We also find that the "age" distribution, 
i.e. the probability distribution of the number of times 
up to time T, that a given site has been a part of the 
shear band (and has hence been changed), has many sim- 
ilarities with models of fragmentation studied in various 
contexts fl7f . 

There has been recently an upsurge of interest in sys- 
tems exhibiting an anomalously slow relaxation. Such 
a behavior is generically reminiscent of a glassy behav- 
ior, and this analogy has motivated a number of stud- 
ies [[l8|. Just to mention one example related to granu- 
lar media, the slow compaction of sand under repeated 
tapping]!^, ^| displays analogies with glasses obtained 
at different cooling rates. A number of different mod- 
eling^ of this compaction process have been proposed 

1 g [§. 



|2g £4|, £5|, ^6|, |27|, [2^, |2£ 

Some of these models emphasize the role of a broad dis- 
tribution of energy barriers which have to be overcome 
through thermal activation. This naturally leads to the 
occurrence of a wide distribution of characteristic times, 
with a slower and slower dynamics as the easiest barriers 
are exhausted. Models of this sort have been looked at 
in a wide variety of contexts ranging from trap models 
p3| and anomalous diffusion in the presence of quenched 
disorder pif to constrained spin system s p5[ , granular 
and ageing in soft solids ]36[. Other ap- 



proaches put more emphasis on the collective nature of 
the necessary rearrangement allowing for a relaxation 
@ H |22l M H |25 



26|, |7|, ||, H U H, £2|. As 

time passes, the relaxation has to become more and more 
cooperative, and hence the barrier is more entropic than 
energetic. Models with entropic barriers have been well 
studied in other contexts too, such as the Backgammon 
Model ||. 

In the model we study in this paper, as we shall see, we 
observe a very slow dynamics which can indeed be com- 
pared to such glassy behavior. We do not include any 
temperature stricto sensu, however, the randomness of 
the local densities can in some way be compared to ther- 
mal noise. The crude classification we proposed above 
between energetic and entropic emphasis is not quite 
suited to our model, where both aspects are simultane- 
ously present. The necessary cooperative nature of effi- 
cient events is included in the search for a minimal path 
where all sites contribute with the same weight. How- 
ever, a local dense configuration can occur at any time, 
and remain quenched thereafter for very long. This is 
like an energy barrier since in order for the minimal path 
to go through this region, all minimal paths with smaller 
energies need to be eliminated This is thus a very rare 
event with the probability becoming smaller and smaller 
as time passes. 

To push forward the analogy with a glassy system, we 
will see that we observe a breakdown of ergodicity, in the 
sense that the activity is not spread uniformly through- 
out the system. Hence if we partition a system into two 
sub parts (even for large sizes), the relative "age" of the 
two subsystems will tend to a broad distribution, and not 
to a narrow one as expected for homogeneous systems. 
This implies that the scaling of the compaction in both 
space and time is expected to be non trivial. 



compaction 



IV. NUMERICAL RESULTS 

In this section we present briefly the most important 
numerical results on both the hierarchical and Euclidean 
lattices. 

The most important quantity of the system is the av- 
erage density. We define it as the mean density of all 
sites not belonging to the minimal path and we denote 
it by (g)(t). The importance of not including the min- 
imal path in the average density is that this definition 
ensures that (g) monotonically increases with time. Fur- 
thermore, as we will see, at late times the minimal path 
mostly remains unchanged. Since we keep refreshing the 
same bonds again and again, the density along the min- 
imal path is simply taken from the known distribution 
Pt{q) and there is no need to incorporate this into (g). 

In our case, as pi — p r is a uniform distribution be- 
tween and 1, it is clear from the rules that the system 
evolves towards the limiting state of g(x, y) = 1 every- 
where. It is natural thus to plot 1 — (g) (t) as done on 

Fig. a 
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FIG. 2: The difference of the average density from its asymp- 
totic value 1 as a function of time. The system sizes are 
L — 32, 64, 128, 256, 512 for the square lattice (upper) and 
N = 2 to N = 7 for the hierarchical lattice (lower) from bot- 
tom to top respectively. The average was done over all the 
inactive sites in the lattice and for an ensemble of 20 to 1000 
samples. System size increases from bottom to top in both 
cases. 



On Fig. H we present both the L x L square (upper) 
and hierarchical (lower) numerical results for (g)(t). We 
can make two immediate observations: The t/L (t/2 N ) 
scaling works nicely up to about unity after which a sys- 
tem size dependent relaxation is observed which is slower 
for larger systems. The density decay seems to be slower 
than any power-law. The other quantity that we study 
in detail is the Hamming distance, i.e., the number of 
different bonds between consecutive minimal paths. We 
denote this quantity by d. The value of d may vary from 
to L (2 N ). 

As can be seen in Fig. ||, in both lattices, at an early 
stage the mean Hamming distance is close to the system 
size (i.e. two consecutive paths do not overlap at all). It 
then decreases monotonically to 0. We recall that when 
the distance is equal to 0, then the two successive con- 
formations of the minimal path are identical, in spite of 
the total renewal of random densities along them. This 
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FIG. 3: The average Hamming distance versus time for square 
(upper) and hierarchical lattice (lower). The same system 
sizes were scaled together as on Fig. ^. On both figures, the 
analytical prediction l/(t + 1) is plotted over the data. Note 
that scaling with system size is excellent for the hierarchical 
lattice while it seems to display systematic corrections for the 
Euclidean case. 



indicates that minimal paths have a tendency to remain 
more and more persistent as the system "ages" . 



V. ANALYTICAL RESULTS ON THE 
HIERARCHICAL DIAMOND LATTICE 

In the following we show how some of the above listed 
properties of our model can be understood analytically 
on the hierarchical lattice. The techniques we use are 
essentially those of extreme value statistics. 



A. Summary of the solution 

We use the hierarchical nature of the diamond lattice 
to calculate the average density for level N knowing the 
results on level N — 1. 
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Let us introduce the following notations: The aver- 
age density of the inactive sites in a (sub)lattice level N 
is Q^ N \t), the density of the active sites is x^ N \t). In 
(sub)lattices we define the active site to be the minimal 
path of this (sub)lattice regardless of whether it is a part 
of the global minimal path or not. We will use the indices 
I, r, u, d for left, right, up, down respectively to indicate 
the parts of a diamond [(sub)lattice] corresponding to the 
illustration on Fig. [j] b,c). 

A generation N lattice is constructed by putting to- 
gether 4 generation N — 1 sublattices (see Fig. with 
two series couplings and one parallel coupling. The series 
coupling of the sublattices is easily taken into account, 
both the density of the minimal path as well as the den- 
sity of the sites in the bulk are simply averaged: 



x (S(N-i )){t) = 1 +4"-i)( t )) = a.("-i)( t ) 

g^ N -^(t) = \ {gf- X \t) + Q^-'Ht)) = g^(t) 

(2) 

where the superscript S refers to two systems coupled in 
series. 

The next step is the coupling in parallel of two series 
couplets. The density of the shear band is simply the 
minimum of that of the two subsystems. The average 
density of the bulk contains two contributions: the av- 
erage densities of the subsystems as well as the average 
density of one of the active paths (the one which is the 
larger of the two contenders for the global minimal path). 



t. Thus we rewrite the second equation of (||) 



,(N) 



(t) = min | : 



.S(N-i) ( , \ S(N- 



- x \u)} 



{Ql^XtJ + gf^iU)} 



(A N - 2 N+1 ) 
2 ( 4 JV _ 2N) 
off , 

(3) 



(4 N - 2 N ) 



where A N is the total number of bonds, 2 N is the num- 
ber of the bonds in a path on a generation N lattice 
and the multiplicative factors in the above equation are 
the appropriate fractions of bonds at generations N (see 
Fig. |). 

There is a further subtlety here. The time t counts the 
total number of updates at generation N. However the 
"time" relevant for a subsystem at generation N — 1 is 
simply the number of times the subsystem itself has been 
updated. Since only one of the two systems in parallel 
is updated at every instant, the "age" of a subsystem at 
level N — 1 is less than t and is denoted by t u and td in 
the above equation. In Appendix H we prove that the 
relative age of either subsystems t u Jt or tdjt is uniformly 
distributed between and 1 in the limit of a large time 



g {N \t) = 



( 4 JV _ 2 iv+i) 1 

{4 N ~2 N ) 1 
1 



(JV-1) 



(t')dt'- 



(4 W - 2 N ) t 



o 

(N- 



1 5(t') J a; (JV - 1) (t-t'))d* / - 



n (4) 
The second term in Eq. |J comes from the competi- 
tion between the minimal paths in the two subsystems 
coupled in parallel. Only one of these is the global mini- 
mum and the larger has hence to be incorporated into the 
density of the system. At every timestep that a subsys- 
tem is updated the minimal path of that subsystem can 
switch to either side. Since the two parallel subsystems 
are entirely disjunct the path changes sides if the mean of 
the random numbers generated along the minimal path 
is larger than the minimal path in the other subsystem. 
This competition is present at all levels of the hierarchy. 

The problem of the minimal path in the parallel cou- 
pling can be thus described by a simple model that we 
call the two site model. The two site model is defined 
as follows: There are two sites, each with a single value 
generated by a random number drawn from a given dis- 
tribution p(x) (in our case at level N it is the sum of 2 N 
independent random numbers each of which is taken from 
the uniform distribution between and 1). We choose the 
site with the smaller value and refresh it with a random 
number generated from the same distribution p8j. The 
dynamics consists of repeating this procedure. Impor- 
tant features of this problem turn out not to depend on 
the distribution p[x), since the entire evolution is only 
based on the ordering of the values. As a result one can 
map any bounded distribution onto a uniform one and 
preserve the same history of the activity. It is thus easy 
to deduce that the probability of having an active site 
in one subsystem for a given time, knowing the age of 
the system, is independent of p. We present in Appendix 
[b| an analytical derivation of relevant properties of this 
problem. 



B. Age distribution 

We have seen that as a result of the parallel coupling, 
the time spent in one subsystem, or the "age" of a sub- 
system, 6, differs from the actual time t and that 0/t 
is uniformly distributed between and 1 in the two-site 
model. Repeating the above argument from the entire 
system down to a single bond, we can extract the statis- 
tical distribution of ages relative to the total time. 

Let pn(T; t) be the statistical distribution that a given 
bond was updated exactly T times at time t in a lattice of 
generation N . Using the above argument, we can relate 
these distributions of different generations through the 
relation: 



p N+ i(T:t) 



PN (t'-t) 



dt' 



(5) 
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with p (T;t) = 6{t - T) and N > 1. In other words 
the age of the bonds or subsystems can be obtained by a 
"fragmentation" process: At every level to get the age of 
the upper and lower arms in the parallel coupling the age 
of the diamond is cut into two pieces with a uniform dis- 
tribution. Not surprisingly equations similar to the above 
are well known in the context of models of fragmentation 



The solution of the above recursion is 



PN(T;t) 



1 (log(t/T)) 



JV-l 



* {N-l)\ 



(6) 



for T < t. Introducing the relative age 9 = T/t, we 
observe that the above distribution becomes independent 
of the time t (the l/t prefactor is absorbed in the measure 
d6 = dT/t). 

We note here that models of fragmentation which are 
described by similar equations usually look for a steady 
state solution i.e. an N independent solution at late 
times. However, in our case, as explained below, the 
N dependence is crucial and has necessarily to be kept. 
Further the order in which N and t are taken to infinity 
is very important as well. 

It is interesting to note that the above distribution can 
be simply expressed in the framework of multifractality, 
which was introduced to characterize the scale depen- 
dence of statistical distributions. This analysis naturally 
provides a generalized "dimensional analysis" of a local 
quantity x, with a distribution Pl{x). We introduce the 
scaling index a and associated fractal dimension f(a) of 
the support of the set of x values defined through 



x oc L a 

xp L (x) cx £/(«)-<* 



(7) 



where d is the space dimension. Alternatively a — 
log(x)/ log(L) and f(a) = d + \og(xp L {x))/ log(L). In 
our case, the local quantity x is the relative age 9 = T/t 
and d = 2, thus 



= log(fl) 
log(L) 

/(a) = 2 + a4 



ii)-io^)){ log( ^ )+ (8) 



+1 + log(log(2)) - 



log(log(L/2)) ] 
log(L) J 



where we have used the Stirling formula, assuming 2 N 
L 1. In this limit, we have 

ft a ) = a | log (- Q ) | l + 21og(2)+log(log(2)) 



-O 



log(2) 

(bi(Z) 



log(2) 




FIG. 4: Multifractal spectrum of the (relative) age distribu- 
tion in the hierarchical lattice, a gives the scaling exponent of 
the age with the system size, and / the corresponding fractal 
dimension of the support of the set of sites contributing to a 
given a. 



where in the limit of an infinite system size, L — > oo, the 
correction term C(l/log(L)) vanishes. Due to this for- 
malism we arrive at a system size independent description 
of the distribution of relative age, although the distribu- 
tion itself depends on L. Moreover, the interpretation of 
the formalism is rather natural. The subset of sites whose 
age scales as a power-law of the system size 8 oc L a has 
a fractal dimension f(a). 

Figure ^ shows the asymptotic form of the multifrac- 
tal spectrum. The range of a values corresponding to a 
positive fractal dimension is a m i n » —5.33 and a max ~ 
—0.15. The scaling exponent characterizing the maxi- 
mum number of sites ao is the one for which / is max- 
imum, i.e. / = 2, and hence ao = — l/log(2) —1.44. 
Let us emphasize that this description is only valid for 
very large times. Otherwise, the finite cut-off in the time 
distribution will affect the multifractal spectrum. More- 
over, we have discarded correction terms which will dis- 
appear as l/log(L), i.e. very slowly. This may render 
this spectrum difficult to observe numerically. 

This analysis shows that the relative age 9 does not 
scale with L in a unique fashion. When computing a mo- 
ment of order m, only one scaling set dominates. The 
precise value of this dominant a depends on m. It corre- 
sponds to the condition d/(a)/da = —m or 



a(m) 



-1 



(m + 1) log(2) 



(10) 



unless the corresponding value of / is negative. The mo- 
ment then varies as 



(m) 



where 



(9) 



l ) ocL r 



t(to) = f(a(m)) — d + ma(m) 



(11) 



(12) 
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In our example, 



r(m) 



log(m + 1) 
log(2) ■ 



(13) 



C. Average density 

Now we can use the above results to get the final form 
of the time evolution of the average density. 

The recursion relation in Eq. is composed of two 
terms. Let us study the first term. As the integration 
operation is additive we can consider separately all the 
components of a lattice of generation iV; from subsystems 
of generation N — 1 right upto individual bonds. 

First we consider the generation 1 lattice, the simple 
diamond. The two site model gives the exact time depen- 
dence of the inactive bonds in the asymptotic limit (see 
Appendix [b]) which is 1 — (t) cx 1 j\[t. In order to get 
the contribution of these bonds we have to complete the 
integral of Eq. |], the expectation value of T -1 / 2 , with 
the correct "age" distribution of these subsystems. Note 
that we calculate the age of a diamond (level 1 object) 
not a bond. Therefore in a level N system we shall use 
PN-i{T]t) for the age distribution. 

1-{qi) <*t~ 1/2 J e l ' 2 p N ^{T-t) dlog(0) 

Jo (N-2)\ &K ' (14) 



= t -l/2 2 N-1 

= r l / 2 L/2 



oo x N-2 
W^Y- 



exp(— x)&x 



where (gj) is the contribution of sublattices of generation 
i to the average density. 

The above result has two important implications. 
First, the time dependence of the sublattice, 1/V~t, is 
preserved on the global scale. Secondly, the statistical 
distribution of ages gives rise to a system size depen- 
dence, i.e. a power-law of L, which in the above example 
displays a trivial exponent 1. More generally, this expo- 
nent is t(— 1/2) as derived above. 

The above expression accounts for about half of the 
bonds. The next term which enters in the coupling is the 
inactive minimal paths in the level 2 subsystems (Fig. [j] 
c). The length of these paths is 2 2 = 4 bonds. The two- 
site model predicts an asymptotic i -1 / 4 time dependence 
for 1 — B( 4 '(t), Thus here we have to use the moment of 
order —1/4 and the age at generation N — 2: 



:i.-(ea) OCT 1 '* I e l / A p N . 2 {T;t) dlog(0) 

oo -JV-3 



,-1/4 



(4/3) 



N-2 



/o (iV-3)! 
= £~ 1/4 (3/4) 2 £ lo g( 4 /3)/l°g(2) 



exp(— x)dx 
(15) 



We see in this example a non-trivial scaling with the 
system size and a slower time dependece. 

We can carry out this same procedure for higher gen- 
ration of subsystems. The length of the path in a level i 
sublattice is n = 2 % thus their contribution is: 

1 - (ft) OC t- l ' n {l - l/ n )*£-l°g(l-l/»)/l°g(2) ( 16 ) 

Thus we observe that the scaling of the mean density 
can be cast into the form of a sum of power-laws with 
a vanishing exponent 1/L = 2~ N and thus a slower and 
slower decay to zero. Each of these terms has a prefac- 
tor which exhibits a different scaling with L, and hence, 
the aging of different system sizes cannot be accounted 
for by a simple reduced time such as t/L. The latter 
only holds for the first subset (half of the system size), 
and not the successive hierarchy of minimal paths. This 
argument explains why the time evolution of the mean 
density seemed to follow a unique curve when plotted 
as a function of t/L for early times. However, as time 
increases, we note a breakdown of this simple scaling, 
and larger systems shows a slower and slower increase 
in the average density. It is interesting to note however 
that for (extremely) large times, (i.e. vanishing exponent 
1/n), the moment will depend only on the combination 
t/L 1 / 10 ^ 2 ). This exponent which appears in L is the a 
value of the largest fractal dimension, / = 2, a = ao in 
the multifractal spectrum. 

The above analysis is, however, valid only for very 
late times, after a long transient. The sum of m iden- 
tically distributed random variables (when each individ- 
ual variable is taken from the uniform distribution) can 
be approximated by a power-law distribution, only at 
late times. Before that, it is well approximated by a 
Gaussian, by the central limit theorem. It is only when 
we are pushed to the tails of the distribution that the 
power-law regime occurs. However the extremising rule 
makes this inevitable, though after a long transient. For 
instance for a generation 3 minimal path, consisting of 
8 bonds in series, this transient ageing period lasts for 
about t ~ 8! ss 40000 time steps, t thus has to be much 
larger than this so that the lower limit of the integral in 
Eq. [l4| can be taken to 0. We see that our computation 
becomes strictly applicable only for extremely late times. 

Finally we put together all the information we have on 
the increase of the density of the inactive bonds of the 
hierarchical lattice i.e.: 

1. The lattice is a collection of two site systems from 
level 1 < i < N where the "ageing" of a bond in 
any level is given by Eq. B7; 



2. The number of bonds in each level gives a prefactor 
of 2 N -y(2 N - 1). 

3. The age distribution at each level results in an ad- 
ditional factor as in Eq. |l6|. 

We finally get the following result for the average den- 
sity: 
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0.01 



io" 1 1 io 1 io 2 io 3 io 4 io 5 



t/2 



N 



along the minimal path in turn plays a crucial role in the 
very slow decay of the mean density in spite of the vanish- 
ing fraction of bonds involved. The departure from 1/2 
varies roughly as y/log(t). Taking such a form into ac- 
count, we see that the average density does not converge 
to 1 any longer, just as if some bonds were quenched close 
to their average value 1/2, up to a very slowly evolving 
correction. Thus numerically, one can achieve a reason- 
able fit of the evolution of q to values different from 1. 
However, as the time window is enlarged, the effective 
asymptotic q increases. Reciprocally, extending the sys- 
tem size, this asymptote decreases. Thus, in spite of the 
quality of the fits which can be produced this way, we un- 
derline the fact that such an approach is only applicable 
to a fixed time or system size window. 



FIG. 5: The test of the analytical result. The numerical data 
are plotted with symbols the corresponding analytical with 
solid lines. The hierarchical level is N = 1, 2, 3, 4 from 
bottom to top respectively. 



,(A0 



(*) = i-£ 



/ 2 N ~ l \ a(i)L z ^ 



(17) 



where (using n = 2 % and L — 2 for a lattice of level TV) 



a(i)= T{l/n)(n\) l ' n n- 2 {l - 1/n) 1 
b(i) = 1/n , > 

log(l - 1/n) ^> 

Z{1) = " log(2) 

Figure |5| visualizes this result compared to the numer- 
ical data. 



D. Intermediate time behavior of the average 
density 

For a minimal path of length I = 2™ , we are interested 
in the maximum value of the sum of i random numbers 
over a number of realizations equal to the age t. For 
large £, the distribution of the average element in the 
sum converges towards a Gaussian of average 1/2 and 
standard deviation 1 / s/Vli. The expectation value of 
the largest such element over a time t is thus such that 



1 



2V31 
V2^ 



exp (-U{x - 1/2) 2 ) dx w - (19) 



It is important to note that this expression is valid for 
large £ and moderate t, whereas we previously considered 
the limit of large t and moderate L. The order of the lim- 
its plays a crucial role. The mean value of the densities 
along the minimal path thus departs only very slowly 
from 1/2. This slow change of the density of the sites 



E. Hamming distance 

Let us now consider the overlap function shown in 
Fig. [| For the hierarchical lattice, the overlap has a sim- 
ple interpretation. We have seen that, at least for large 
times, most of the activity essentially takes place along 
the same path. However, from time to time, the minimal 
path jumps from one conformation to another, whose dis- 
tance to the previous one is quantified by the Hamming 
distance. The scarcity of the jumps allows us to neglect 
the occurrence of simultaneous multiple jumps. Let us 
define the probability, Pi,{d,t), as the probability that a 
jump equal to d = 2" takes place, i.e. the probability 
that the current path differs from the previous one by d 
sites in a system of size L = 2 N . This quantity, apart 
from containing information about how the average value 
of d changes with time, is also the natural analogue of an 
"avalanche distribution" in this model. As will be seen 
further down, this quantity does indeed decay for large 
times as a power-law of the distance d like in many other 
self-organized critical models. However the distribution 
has a time-dependent prefactor unlike other models with 
a true steady state. 

For n = N, we have to consider a jump at the largest 
scale available in the system. At this level, the lattice can 
be coarse-grained as a generation 1 lattice. The proba- 
bility for such a jump to occur is equal to the probability 
that in a two-site model, the activity moves from one 
site to the other one. We show in the appendix that this 
probability, is equal to l/{t+ 1), and thus for large times, 



P L (d = 2 J \t) a 1/t 



(20) 



Let us now consider a smaller jump size, i.e. n = 
N — 1 . This means that one half of the actual minimal 
path should move to a different configuration. Thus we 
focus on a subsystem of size L/2, whose age is T. In the 
appendix, we show that the probability for such a n ag e 
is q(T,t) = 2(T+l)/(t+l)(i + 2) « 2T/t 2 , (see Eq.|Af). 
Moreover, we have two such subsystems in series and thus 
the probability that the Hamming distance is L/2 in a 
system of size L is approximately twice the probability 
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d 

FIG. 6: P(d, t)t scaled by the system size for the hierarchical 
lattice of generation N = 6 at times t/L = 100, 200, 300 and 
1000 and for generation N = 7 for times t/L = 100 and 500. 
Time is measured in terms of the system size and the straight 
line has a slope —2. 

that in one subsystem the Hamming distance is equal to 
the system size, P^/ 2 (L/2,t). Integrating over all times 
T with the above probability we have 

r* 2T 1 

P L (L/2,t)^2j o _-dT = 4/t (21) 

For smaller jumps, we can repeat the same argument 
recursively, to obtain 



x N-n+l 
2 2(N-n) L 2 
t = Id? 

As one can see this result agrees well with the numer- 
ical results shown on Fig. ^. Moreover, this expression is 
to be compared with the Hamming distribution obtained 
with a logarithmic measure for d on the Euclidean lat- 
tice [||, which has the same functional form. Hence this 
distribution is the same for the two lattices despite their 
connectivities being very different. Note also that this 
quantity scales perfectly with the system size and is de- 
scribed at all times by the reduced time distribution t/L 



as evidenced in Fig. |3| Thus the behaviour of the Ham- 
ming distance is much simpler than the slow density in- 
crease in the system. 

VI. DISCUSSION AND CONCLUSION 

We have presented simulation results and an asymp- 
totic analysis of the behavior of the the optimization and 
restructuring model on the hierarchical lattice. The two 
lattices are very different in structure, yet they also ex- 
hibit remarkably similar features. 

The Hierarchical lattice is easier to analyse due to its 
recursive structure. For example, the very specific con- 
nectivity of sites on this lattice ensures that large jumps 
are always possible, though rare at late times. In the 
Euclidean lattice, these are strongly suppressed by a fur- 
ther feedback effect: the localization of the path limits 
the density increasing effect of the dynamics to a small 
region around the path which in turn intensifies the lo- 
calization. This feature also results in the density map 
being very different in the two cases. In the Euclidean 
case [|] the inhomogeneities in a late-time snap shot of 
the system are much more enhanced. Another difference 
is that changes in conformation in the hierarchical lat- 
tice are organized in a strictly hierarchical way. This is 
clearly not so on the square lattice, where randomness 
and self-organization play an important role. Neverthe- 
less, the overall behavior of the two type of lattices is 
remarkably similar. 

These similarities are most apparent in the time and 
size dependence of the average density (Fig. ||). In both 
cases we have a data collapse for short times, while for 
longer times the dynamics becomes slower and slower as 
the size of the system increases. For the hierarchical lat- 
tice, we have obtained an analytic expression for the scal- 
ing of the Hamming distance, for the local age distribu- 
tion, and its multifractal spectrum, and the asymptotic 
average density evolution with time. In particular, the 
mechanism behind the breakdown of ergodicity, and the 
unusual size dependence of the density evolution can be 
traced back to the multifractal distribution of age. The 
latter provides a novel scenario for "glassy" ageing. 
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APPENDIX A: TIME PARTITIONING IN THE 
TWO-SITE PROBLEM 

We consider here the two-site model and prove that 
after the elapse of a time t, the probability that a given 
site has been visited T times is uniformly distributed 
between and t. The following proof is valid for any 
distribution p(x). Let us compute q(T, t), the probability 
that the active site has been refreshed T times up to time 
t. Let the value of the recently refreshed site be denoted 
by x and the inactive site at t be bt . At time t two things 
may happen: 



i) either x < &t_i, and thus T t = 
happens with probability 
(* + !)/(* + 2) 



1. This 



ii) or x > bt- 
l/(t + 2). 



and thus T t — t—Tt—i with probability 



Note that for any distribution p(x) the activity change 
can only be due to the fact that the largest generated 
random number up to time t is at instant t. This happens 
with probability l/(t+ 2) because in an independent time 
series of random numbers the largest number is equally 
likely to be anywhere. At time t = we have to initialize 
the system by generating two random numbers for the 
two sites. This is the reason for the shift in time from t 
tot + 2. 

Now we can write down a simple evolution equation 
for q(T,t) 

q(T,t) = g(T-l,*-l)^+g(i-T-M-l)^ (Al) 
The general solution of the recursion can be written as 



q(T,t) = A 



(T + l) 



(i + l)(t + 2) 



B 



(A2) 



It is simple to compute q(t, t) from the above recurrence, 
and get 



q{t,t) 



(t + 2) 



(A3) 
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and thus B = and A = 2. Thus finally 

2(T + 1) 



q(T,t) 



(t + l)(t + 2) 



(A4) 



Now, the number of updates of the other site is simply 
q(t — T, i), thus we can formulate the probability distri- 
bution that a site has been updated T times: 



r(T,t) = ^(q(T,t)+q(t-T,t)) = 
(T+l) + (t-T + l) 



(A5) 



+ l)(t + 2) (t + 1) 

which is independent of T. Thus r(T,t) is uniform 



APPENDIX B: AVERAGE DENSITY OF THE 
INACTIVE SITE IN THE TWO SITE MODEL 




100 1000 10000 

t 



FIG. 7: The test of the results in Eq. B7 against numerical 



simulations on the same model. The solid lines are the an- 
alytical solutions, the symbols indicate numerical simulation 
results. The system sizes are n = 2(o), 4(+), 8(D), 16(x) 



Since the hierarchical lattice can be considered as a set 
of two-site systems, we need to consider the distribution 
of the maximum in a two site model in the case when each 
of the sites is taken from a distribution p n (x). Here the 
subscript n denotes that this is the distribution for a sum 
of n independent random numbers each of which is taken 
from the uniform distribution. We need only consider 
the case when n — 2 N as the length of the hierarchical 
lattice can only be of this form. 

Unfortunately, p n (x) is difficult to formulate in a gen- 
eral way. For large n this is a Gaussian for moderate val- 
ues of x. However very close to the extremes and 1, it is 
a power-law as we will see below. It is this regime which 
is asymptotically reached and hence relevant for our pur- 
poses. We hence only consider the regimes x < 1/n or 
x > 1 — 1/n (cases when only one number out of n may 
reach its extreme value 1 in the x < 1/n case and in 
the x > 1 — 1/n case). 

Let us recall the formula for the average value of the 
largest generated number up to time t when each of the 
individual numbers x is taken from a distribution p n {x): 

B (n \t) = (t + 1) f x Pn (x)P^(x)dx, (Bl) 
Jo 

where P n {x) is the cumulative distribution of p n (x). 
Thus Pn(x) accounts for the probability that the t other 
numbers are less than x. The t + 1 factor is needed to 
take into account the fact that the position of the largest 
number can be anywhere in time. The index n indicates 
that the distributions describe the average of n indepen- 
dent, uniformly distributed random numbers. 

For ( > 1 we have Pn(x) <C 1 for most values of x, 
except for a 1/n neighborhood of 1. 

So in the integration the most important contribution 
comes from the part that is close to 1. This permits us 
to restrict the integral to the part that we can calculate 



without loss of consistency: 

B^(t) ~ f (t + ^xp^P^dx (t » 1) 

JX-l/n 

(B2) 

The probability distribution close to the limits takes the 
following forms 

n n 

Pn{x) x<1 / n = -TT X"- 1 (B3) 

(n — 1)1 

n n 

Pn(x)\ x> i-y n - J^Jyi 1 -^ 1 ( B4 ) 
The cumulative distribution is the integral of the above: 

n n-l 

P n (x)\x>i-i/ n = 1 - (n _ jT; (1 - X T ( B5 ) 

Let us now turn back to Eq. |B^. Using an x = 1 — y 
variable replacement and doing integration in parts we 
arrive at the following formula after neglecting the expo- 
nentially decaying parts: 

B {n \t) = (t + 1)- / (l-y)y n - l x 

fa - 1)! Jo 

n'-fa^J dy= 

/•1/n / n-1 \ *+! 

We rewrite the integrand in a (• • -) t+1 = exp{(i + 
1) log(- ■ ■)} form and make a Taylor expansion in y n 
around y — to the second order. The result can be 
written in the following form that we use in our calcula- 
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tions: On Fig. |?] we can see that the above approximation is 

excellent for small n and large t. 



B^(t) = 1 - r(l/n)(n!) 1 /"n- 2 ft + 

+U { t l/n+l I 



-1/n 



(B7) 



